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Foreword 


Hydrographic  and  bathymetric  charts  of  the  world’s  coastal  waters 
are  essential  to  the  missions  of  many  U.S.  Navy  programs.  A  coastal 
survey  requirements  study  by  the  Defense  Mapping  Agency  (DMA) 
showed  a  200-  to  300-year  backlog  if  the  necessary  bathymetry 
data  are  to  be  gathered  by  conventional  ship  surveying  methods. 

The  Mapping,  Charting,  and  Geodesy  (MC&G)  Division  of  the 
Naval  Oceanographic  and  Atmospheric  Research  Laboratory 
(NOARL)  was  tasked  to  develop  the  Airborne  Bathymetric 
Survey  (ABS)  system  and  its  supporting  software.  In  response,  the 
MC&G  Division  produced  the  ABS  system  design,  the  NOARL 
Scanner,  and  a  laser  sounder.  These  instruments,  together  with  a 
global  positioning  system  and  an  inertial  navigation  system  aboard 
a  U.S.  Navy  P-3  aircraft,  are  intended  to  survey  clear,  shallow 
coastal  waters.  In  addition,  software  was  written  to  process  the 
multispectral  data  provided  by  the  NOARL  Scanner  and  to  calculate 
the  shallow-water  bathymetry  based  on  laser-  and  boat-derived 
soundings. 

The  instrumentation,  the  software,  and  the  algorithms  that  resulted 
from  this  effort  are  described  in  this  report.  This  advanced  technology 
could  reduce  the  DMA  survey  data  backlog  and  help  provide  current, 
high-resolution  hydrographic  and  bathymetric  data  to  users. 

i— 

W.  B.  Moseley  R.  Elliott,  Commander,  USN 

Technical  Director  Commanding  Officer 
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Executive  Summary 


The  Mapping,  Charting,  and  Geodesy  (MC&G)  Division  of  the  Naval 
Oceanographic  and  Atmospheric  Research  Laboratory’s  Ocean  Science 
Directorate  is  the  primary  activity  within  the  U.S.  Navy  for  conducting  research 
and  development  in  direct  support  of  naval  MC&G  requirements.  The  Mapping 
Sciences  Branch  of  the  MC&G  Division  was  tasked  to  develop  the  algorithms 
and  software  necessary  to  process  data  and  to  produce  bathymetry  from  the 
multispectral  scanner  of  the  Airborne  Bathymetric  Survey  system  when  it 
becomes  fully  functional.  The  software  developed  for  this  system  is  called  the 
Multispectral  Image  Depth  Analysis  System,  or  MIDAS.  The  software 
processes  multispectral  data  in  conjunction  with  laser-  or  boat-derived  soundings 
to  produce  high-resolution  bathymetric  grids.  The  system  is  designed  to 
estimate  depths  in  clear,  shallow  coastal  waters  down  to  20  m.  The  depth' 
accuracies,  derived  from  the  system  using  boat  soundings  for  control,  ranged 
from  0.3  m  to  1.4  m  in  waters  extending  to  3  m  and  10  m,  respectively.  The 
horizontal  positioning  accuracy  of  the  system  was  in  the  25 -m  range,  as 
determined  from  a  ground  survey. 
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Multispectral  Software  Development  for  the 
Airborne  Bathymetric  Survey  System 


Introduction 

This  report  presents  the  results  of  the  scanner 
software  and  algorithm  development  for  the  Airborne 
Bathymetric  Survey  (ABS)  System  developed  by 
the  Naval  Oceanographic  and  Atmospheric  Research 
Laboratory  (NOARL).*  A  brief  background  of  the 
development  of  the  system  and  its  scope  and 
objectives  is  given,  followed  by  the  results  of  the 
scanner  software  effort. 

The  ABS  system  comprises  a  hydrographic 
airborne  laser  sounder  (HALS),  a  multispectral 
scanner  (MSS),  a  global  positioning  system  (GPS), 
and  an  inertial  navigation  system  (INS)  on  board  a 
Navy  P-3  hydrographic  airplane.  The  Naval  Air 
Development  Center  (NADC)  installed  the  hardware 
aboard  the  P-3.  Figure  1  depicts  the  positioning  of 
the  system  components  on  board  the  aircraft.  The 
NOARL  Scanner  (NS)  is  located  on  the  center  line 
of  the  plane  aft  of  the  laser  sounder  system.  The 
system  is  designed  to  collect  hydrographic  infor¬ 
mation  over  clear,  coastal  waters  at  depths  up  to 
20  m. 

The  ABS  system  resulted  from  three  develop¬ 
ments:  a  requirement  for  hydrographic  and 
bathymetric  information  from  the  world’s  coastal 
areas  that  is  essential  to  the  mission  of  many  Navy 
programs;  a  Defense  Mapping  Agency  (1978) 
coastal  survey  requirements  study  that  showed  a  con-; 
ventional  survey-ship  backlog  of  200  to  300  years, 
or  8250  nmi2;  and  advances  in  laser  and  satellite 
technology  that  made  the  ABS  system  an  achievable 
goal  (Polcyn,  1976;  Enabnit,  1979;  Compton,  1988). 

The  system  was  developed  in  two  phases.  Phase  I 
MSS  development  was  to  assemble  and  test  a  nine- 
channel  scanner  and  data  processing  algorithms  to 
meet  operational  requirements  (discussed  later).  A 
description  of  Phase  I  for  the  multispectral  software 


•Formerly  the  N«v»I  Ocean  Research  and  Development  Activity 
(NORDA). 


is  given  by  Kalcic  et  al.  (1989).  Phase  II  develop¬ 
ment  was  to  enhance  algorithms  and  computer 
software  necessary  to  fully  position  and  process 
scanner  data  to  approach  International  Hydrographic 
Office  (IHO)  standards. 


Operational  Requirements 

ABS  system  accuracy  requirements  must  meet 
IHO  standards  (Table  1)  and  ensure  the  detection 
of  all  hazards.  To  provide  a  high  likelihood  of  hazard 
detection,  the  pixel  size  of  the  MSS  was  required 
to  be  1  m  or  less,  so  that  such  a  feature  would  be 
spread  over  at  least  five  pixels.  The  ABS  system 
was  also  required  to  provide  200%  coverage  of  the 
bottom  because  of  the  possibility  that  noise  would 
mask  subsurface  hazards.  The  Naval  Oceanographic 
Office  (NAVOCEANO)  plans  to  obtain  that  cov¬ 
erage  by  flying  the  system  in  two  successive  flight 
lines  spaced  250  m  apart.  The  line  spacing  and  the 
overlap  require  a  swath  width  of  at  least  500  m. 
This  width  is  achieved  by  having  a  field  of  view  of 
54  degrees  or  larger.  In  addition  to  the  parallel  flight 
lines,  cross-check  lines  will  be  flown  approximately 
2500  m  apart. 


Processing  Requirements 

The  ABS  system  is  expected  to  operate  an 
average  of  160  to  200  hours  per  year.  The  process¬ 
ing  requirements  for  the  NS  data  stem  from  the 
collection  rate  of  100  scanlines  per  second.  Data 
collected  average  36,000  scanlines  per  hour,  or 
7,200,000  scanlines  per  year  or  331.776  Gbytes 
per  year.  To  process  all  die  data,  the  daily  proc¬ 
essing  requires  that  28,800  scanlines,  over  fifty 
512x512  byte  images,  each  with  five  channels, 
be  processed  each  day  (based  on  250  working 
days/year). 
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-34.5’ 

Center  of  gravity 

+  621' 

0 
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Figure  1.  ABS  system  configuration  on  board  a  P-3  aircraft.  The 
displacements  of  the  system  components  from  a  point  in  the  center  of 
the  aircraft  are  given  as  x,  y.  and  z  in  the  table.  The  system  also  has 
recorders  on  board. 


The  NOARL  Scanner 

The  NS  is  a  9-channel  MSS  with  increased 
sensitivity  in  the  blue-green  portion  of  the  spectrum. 
Schematic  representations  and  specifications  for  the 
NS  are  given  in  Appendix  A. 

The  spectral  resolution  of  the  NS  is  given  in 
Table  2.  Bands  1-6  are  visible  light  bands  and  bands 
7-9  are  infrared  bands.  The  NS  spectral  bandwidths 
are  similar  to  National  Aeronautics  and  Space 


Administration’s  (NASA)  Landsat  Thematic 
Mapper  (TM),  except  that  the  first  four  bands  of 
the  NS,  the  water-penetrating  bands,  make  up  the 
first  two  bands  of  the  TM.  NS  bands  5-9  are 
equivalent  to  TM  bands  4-7. 

The  data  are  recorded  on  a  14-track,  high-density 
tape  on  board  the  aircraft.  The  high-density  tapes 
are  decommutated  to  6250-bpi  density  computer 
compatible  tapes  (CCT).  Table  3  shows  the  format 
of  the  data  records  on  the  6250-bpi  tapes. 
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Table  1.  A8S  system  operational  requirements. 


Threshold 

Requirement 

(Based  on  IHO  Standards) 

(Based  on  NAVOCEANO  Requirements) 

Positioning 

15-m  platform  relative  to  shore 

5-m  sensor  relative  to  platform 

Depth  Determination 

±  0.3  m  in  depths  0-30  m 
±  1.0  m  in  depths  30-100  m 
±  1%  in  depths  greater  than  100  m 

Hazard  Detection 

Minimum  Target  Size  is  5-m  diameter 
in  depths  less  than  20  m 

Least  Depth  Determination 

100%  spatial  coverage 

Table  2.  Spectral  bandwidths 
of  the  NOARL  scanner. 


Band 

Spectral  Band 

1 

450  -  480  nm* 

2 

450- 520  nm 

3 

520- 550  nm 

4 

550-600  nm 

5 

630  -  690  nm 

6 

760  -  900  nm 

7 

1 .55  -  1 .75  pm** 

a 

2.08 -2.35  pm 

0 

10.4 -12.5  pm 

*nm  nanometers 
**pm  micrometers 


Multispectral  Image  Depth 
Analysis  System 

The  Multispectral  Image  Depth  Analysis 
System  (MIDAS)  is  a  software  system  that 
processes  the  multispectral  data  to  produce  high- 
resolution  bathymetric  grids.  The  system  consists 
of  1 1 .000  lines  of  Fortran  77  source  code,  making 
up  109  programs  and  subroutines.  Appendix  B 
describes  the  algorithms  and  techniques  that  are 
used  in  MIDAS  software.  Kalcic  et  al.  (1988)  list 
the  source  codes  and  provide  user  documentation 
for  the  software.  These  methods  result  from  much 
of  the  work  done  in  Phase  I  of  the  ABS  system 
development  (Kalcic  et  al.,  1989)  and  continued  in 
Phase  II. 

The  Phase  II  software  implementation  of  the 
algorithms  is  described  in  this  section.  The  soft¬ 
ware  is  divided  into  two  major  processing  areas. 
The  first  process,  the  deblocking  procedure,  divides 
continuous  flight  lines  into  subareas  and  pieprocesses 


the  navigation.  The  second  process  rectifies  the 
multispectral  image  for  attitude  distortion;  regis¬ 
ters  it  to  a  1-m,  equal-area  projection  grid;  maps 
laser-derive  derived  depths  to  the  multispectral 
image;  and  calculates  depths  and  outputs  them  to 
an  image  file.  Figure  2  is  an  operational  diagram 
for  MIDAS. 

Preprocessing 

The  deblocking  procedure,  DEBLOCK,  reads  the 
decommutated  tape,  and  divides  one  flight  line  into 
manageable  subareas.  This  procedure  operates  in 
batch  mode.  The  size  of  the  subareas  is  adjustable, 
but  512  x  624  pixels  or  0.5  x  0.6  km,  was  convenient 
because  most  image  display  devices  are  capable  of 
512  x  512  pixel  display.  Navigation  and  attitude  data 
are  validated  and  a  position  is  computed  for  each 
scanline  by  module  GEONAV.  At  completion,  a 
file  of  subareas  is  generated  for  input  to  the  MIDAS 
processing  system. 

MIDAS  Processing  System 

MIDAS  performs  the  geometric  corrections 
(GEORECT),  maps  the  laser-derived  depths  to 
the  multispectral  image  (GEOMATCH),  performs  the 
depth  calculations  (least-mean-square  (LMS)  or 
DEEP),  and  outputs  a  depth  file.  MIDAS  allows 
the  analyst  to  view  the  data  before  processing  and 
to  view  the  results  after  each  process.  MIDAS  has 
the  option  to  run  in  two  modes:  interactive  or  pro¬ 
duction.  A  flow  diagram  in  shown  in  Figure  3. 

Modes  of  Operation 

Interactive  mode.  This  option  allows  the  analyst 
to  view  the  MSS  data  prior  to  processing,  or  the 
results  at  any  step.  The  analyst  has  the  option  of 
processing  any  subarea  within  the  flight  line  and 


Table  3.  NS  data  format  (1  record  of  6250  bpi  tape). 


Word 

Byte 

Parameter 

Bits 

Description 

1 

0 

Spare 

0-7 

1 

Scan  Line  Count 

0-3 

10,000's 

1 

4-7 

100,000'S 

2 

0-3 

100’s 

2 

-  * 

4-7 

1 ,000’S 

3 

0-3 

1'S 

3 

4-7 

10'S 

2 

0 

Spare 

0-7 

1 

Scan  Speed 

0 

100's  of  scans/sec 

1 

■  ■ 

1-7 

unused 

2 

•  a 

0-3 

1's  of  scans/sec 

2 

e  • 

4-7 

10's  of  scans/sec 

3 

Spare 

0-7 

3 

0 

Roll 

0-3 

6  MSB's’  of  14 

0 

■  m 

6-7 

bit  6  is  unused  validity  bits 

1 

m  m 

0-7 

8  LSB's1  of  14 

2 

Pitch 

0-5 

6  MSB's  of  14 

2 

.  a 

6-7 

bit  6  is  unused  validity  bits 

3 

.  e 

0-7 

8  LSB's  of  14 

4 

0 

Shot  Time 

0-1 

2  MSB's  of  10  of  msec 

0 

m  m 

2-5 

4  bits  of  sec 

0 

m  m 

6-7 

bit  6  is  unused  validity  bits 

1 

•  • 

0-7 

8  LSB's  of  20  of  msec 

2 

Spare 

0-7 

3 

Time 

0-1 

2  MSB's  of  10  of  msec 

3 

•  • 

2-6 

unused 

3 

•  ■ 

7 

MSB  of  17  of  sec 

5 

0 

a  a 

0-7 

8  LSB's  of  10  of  msec 

1 

a  a 

0-7 

bits  6-15  of  17  of  sec 

2 

a  a 

0-7 

8  LSB's  of  17  of  sec 

3 

Altitude 

0-7 

8  MSB's  of  16 

6 

0 

a  a 

0-7 

8  LSB'sF  of  16 

1 

Drift  Angle 

0-7 

8  MSB's  of  16 

2 

a  a 

0-7 

8  LSB’s  of  16 

3 

True  Heading 

0-7 

8  MSB's  of  16 

7 

0 

a  a 

0-7 

8  LSB's  of  16 

1 

Latitude 

0-7 

8  bit  exponent  of  48  bit  FP* 

2 

Latitude 

0-7 

bits  32-39  of  40  bit  mantissa 

3 

a  a 

0-7 

bits  40-47  of  40  bit  mantissa 

8 

0 

a  a 

0  7 

bits  16-23  of  40  bit  mantissa 

1 

a  a 

0-7 

bits  24-31  of  40  bit  mantissa 

2 

a  a 

0-7 

bits  0-6  of  40  bit  mantissa 

3 

Longitude 

0-7 

8  bit  exponent  of  48  bit  FP 

g 

0 

a  a 

0-7 

bits  32039  of  40  bit  mantissa 

1 

a  a 

0-7 

bits  40-47  of  40  bit  mantissa 

2 

a  a 

0-7 

bits  16-23  of  40  bit  mantissa 

3 

a  a 

0-7 

bits  24-31  of  40  bit  mantissa 

10 

0 

a  a 

0-7 

bits  0-6  of  40  bit  mantissa 

1 

Horizontal  Velocity 

0-7 

8  bit  exponent  of  32  bit  FP 

2 

a  a 

0-7 

bits  16-23  of  40  bit  mantissa 

3 

a  a 

0-7 

bits  24-31  of  40  bit  mantissa 

11 

0 

a  a 

0-7 

bits  0-6  of  24  bit  mantissa 

1 

Direction 

0-7 

8  bit  exponent  of  32  bit  FP 

2 

a  a 

0-7 

bits  16-23  of  24  bit  mantissa 

3 

a  a 

0-7 

bits  24-31  of  24  bit  mantissa 

’most  significant  bit  'least  significant  bit  ’floating  point 
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Table  3.  Continued. 


Word 

Byte 

Parameter 

Bits 

Description 

12 

0 

■  ■ 

0-7 

bits  0-6  of  24  bit  mantissa 

1 

GPS  Time  at  Solution 

0-7 

8  bit  exponent  of  48  bit  FP 

2 

•  ■ 

0-7 

bits  32-39  of  40  bit  mantissa 

3 

•  • 

0-7 

bits  40-47  of  40  bit  mantissa 

13 

0 

0-7 

bits  16-23  of  40  bit  mantissa 

1 

■  * 

0-7 

bits  24-31  of  40  bit  mantissa 

2 

•  ■ 

0-7 

bits  0-6  of  40  bit  mantissa 

3 

BB1  Temp 

0-2 

10’s  of  degrees 

3 

■  " 

3 

sign  (0-+) 

3 

4-7 

unused 

14 

0 

■  • 

0-3 

0.1 's  of  degrees 

0 

•  ■ 

4-7 

1’s  of  degrees 

1 

BB2  Temp 

0-2 

10  s  of  degrees 

1 

■  ■ 

3 

sign  (0-+) 

1 

■  ■ 

4-7 

unused 

2 

•  • 

0-3 

0.1 ’s  of  degrees 

2 

■  ■ 

4-7 

I  s  of  degrees 

3 

Image  Data  (700  bytes) 

1  byte  pixels 

IS 

0 

189 

3 

BB1  Radiant  Temperature 

0-7 

190 

0 

Sphere 

0-7 

1 

BB2  Radiant  Temperature 

0-7 

2 

Gain 

0-3 

2 

Channel  ID 

4-7 

3 

Fill  Pattern 

0-7 

55  (hex) 

of  setting  processing  limits  within  the  subarea. 
Setting  processing  limits  within  the  subarea  allows 
the  analyst  to  eliminate  image  portions  that  are 
contaminated  with  sun  glint/glare  or  are  of  poor 
quality  for  other  reasons. 

Production  mode.  This  option  allows  the  analyst 
to  set  initial  parameters — i.e.,  region  of  interest, 
starting  subarea.  Successive  subareas  are  processed 
until  the  end  of  the  flight  line  is  reached. 

Data  Storage 

MIDAS  contains  three  common  data  areas:  static, 
dynamic,  and  scratch.  The  static  area  contains  data 
structures  that  remain  constant  throughout  processing. 
The  dynamic  area  contains  data  structures  that  may 
change  between  modules.  The  scratch  area  is  for 
mapping  local  variables  into  memory.  This  area  of 
memory  is  used  repeatedly  by  various  routines,  so 
it  is  referred  to  as  “scratch”  memory.  All  data, 
including  image  data  for  the  subarea  being  proc¬ 
essed,  reside  in  memory  at  all  times.  This  data 
storage  capability  requires  a  large  program  size 
(e  g.,  12  Mbytes  for  6  multispectral  channels).  The 
tradeoffs  are  that  data  input  and  output  (I/O)  are 


reduced  and  the  storage  of  intermediate  files  need 
not  be  eliminated.  The  display  portion  of  MIDAS 
reads  the  data  directly  from  memory,  which  also 
reduces  I/O. 

Modules  Description 

PNAV.  Module  PNAV  is  part  of  the  deblocking 
process.  PNAV  reads  the  magnetic  tape  file  and 
loads  the  image  and  navigation  data  into  a  common 
memory  area. 

Each  time  PNAV  is  invoked,  a  subarea  is  loaded 
into  the  common  area.  For  the  first  subarea,  scanlines 
are  searched  until  the  millisecond  portion  of  the 
shot-time  is  zero.  This  event  is  defined  as  a  shot¬ 
time  rollover.  This  scanline  contains  the  first 
navigation  position  so  this  is  where  processing  starts. 

As  each  scanline  is  read  the  image  data  are  read 
into  common  memory  in  2-byte  or  16-bit  variables. 
The  extra  space  is  allotted  to  handle  the  Airborne 
Multispectral  Pushbroom  Scanner  (AMPS)  12-bit 
data.  The  navigation  data  unpacking  is  done  at  this 
point.  Only  those  fields  used  in  processing  are 
unpacked.  They  include  latitude,  longitude,  GPS 
heading,  INS  heading,  GPS  time,  shot  ume,  roll, 
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Figure  2.  Operational  diagram  for  MIDAS. 


pitch,  and  drift.  The  image  and  navigation  data  are 
read  until  enough  lines  for  a  subarea  plus  an  overlap 
area  or  the  end  of  the  flight  line  is  encountered. 
The  navigation  data  are  read  until  shot-time  rollover 
is  found.  The  overlap  data  are  stored  in  local  buffers, 
which  are  read  in  when  the  successive  area  is 
processed. 

GEONAV.  Module  GEONAV  is  part  of  the 
deblocking  process.  Aircraft  navigation  and  attitude 
data  are  checked  for  validity.  The  GPS  parameters 
of  latitude,  longitude,  and  heading  are  shifted  back 
2  seconds  to  account  for  the  time  delay  in  receiving 
them.  A  time  and  a  position  are  computed  for  each 
scanline. 


A  data  validity  mask  is  built  by  scanning  the 
GPS  times.  A  GPS  time  of  zero  denotes  a  data 
dropout.  In  the  event  of  a  dropout,  it  is  assumed 
that  all  parameters  received  over  the  slow  bus 
(Kalcic  et  al.,  1989),  which  includes  latitude,  lon¬ 
gitude,  altitude,  drift,  INS  heading,  and  GPS  heading, 
are  invalid.  The  previous  GPS  times  are  scanned 
until  a  valid  GPS  time  is  encountered.  A  valid  time 
is  computed  by  linear  interpolation  as  a  function  of 
scan  rate. 

Altitude,  drift,  and  INS  heading  are  replaced  by 
the  closest  valid  value.  Latitude,  longitude,  and  GPS 
heading  are  valid  only  when  the  millisecond  portion 
of  the  shot  time  goes  to  zero.  The  shot  time  is  a 
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time-keeping  parameter  used  to  flag  GPS  readings. 
The  shot-time  clock  resets  to  zero  whenever  a  GPS 
value  becomes  valid.  The  parameters  are  replaced 
only  when  a  data  dropout  occurs  during  a  rollover. 

Aircraft  attitude  is  scanned  for  noise  spikes.  A 
difference  threshold  of  0.5  degree  for  roll,  pitch 
and  drift  and  a  difference  threshold  of  100  m  for 
altitude  are  used  to  flag  spikes.  When  a  spike  occurs 
the  parameter  is  replaced  by  the  previous  valid  value. 

The  shot  time  received  over  the  fast  bus  (Kalcic 
et  al.,  1989)  is  updated  every  scanline.  If  an  update 
does  not  occur,  then  the  time  is  replaced  by  linear 
interpolation. 

A  time  is  generated  for  each  scanline  using  the 
GPS  time  and  the  shot  time.  The  initial  time  is 
the  GPS  time  at  the  first  shot-time  rollover  of  the 
flight  line.  Successive  times  are  incremented  from 
the  initial  time. 

The  GPS  parameters,  latitude,  longitude,  and 
heading,  are  updated  approximately  every  second. 


With  a  scan  rate  of  100  lines  per  second,  a  positional 
update  occurs  every  100  scan  lines.  These  parameters 
become  valid  at  the  shot-time  rollover.  If  the  first 
block  within  the  flightline  is  being  processed,  then 
it  is  assumed  that  the  first  line  has  a  shot  time  of 
zero.  If  the  first  block  is  not  being  processed,  the 
first  scanline  has  a  valid  position  computed  from 
the  previous  block.  The  data  are  then  scanned  for 
a  second  shot-time  rollover.  The  data  validity  flag 
is  checked.  If  the  data  are  invalid,  then  the  next 
rollover  is  found  with  valid  data.  Up  to  5  seconds 
of  data  can  be  interpolated.  When  two  valid  posi¬ 
tions  are  found,  a  linear  interpolation  is  performed. 

Roll,  pitch,  and  drift  are  updated  every  100  msecs 
or  every  10  scanlines.  These  parameters  are  low 
pass  filtered  using  a  13-point  Martin  filter  with  a 
frequency  cutoff  of  0.05  cycles/0.01  secoiJs. 

GEORECT.  Module  GEORECT  geometrically 
corrects  the  multispectral  image  for  aircraft  attitude 
(roll,  pitch,  and  drift).  Mapping  coefficients  are 
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computed  by  using  the  image  boundaries  and  the 
interior  positions  as  control  points  (see  App.  B). 
Geometric  corrections  are  performed  on  each  pixel. 
The  corrected  pixels  are  mapped  to  a  1-m  grid  via 
the  mapping  equations.  Oversampling  and  changes 
in  the  roll,  pitch,  and  drift  angles  may  cause  more 
than  one  pixel  to  occupy  the  same  cell.  In  these 
cases  the  digital  count  values  are  averaged.  At  this 
point,  land  is  identified  and  mapped  to  the  grid  and 
invalid  data  are  rejected. 

The  distance  between  the  center  of  pixels  increases 
as  the  distance  from  nadir  increases;  i.e.,  at  an 
altitude  of  500  m  this  distance  becomes  1.7  m  at 
the  outermost  edge.  When  sampled  to  a  1-m  grid, 
this  effect  causes  data  gaps  along  the  scanline.  Data 
gaps  between  scan  lines  may  occur,  depending  on 
sampling  rate,  air  speed,  and  changes  in  aircraft 
attitude.  Nearest-neighbor  resampling  is  used  to  fill 
in  these  data  gaps.  Only  those  pixels  that  have 
adjacent  neighbors  are  resampled;  otherwise,  a  depth 
is  not  computed  for  that  position.  The  order  of  the 
nearest  neighbor  search  is  as  follows:  left,  right, 
below,  above,  lower  left,  lower  right,  upper  left, 
upper  right. 

Band  6  (760-900  nm)  is  used  as  a  land  mask 
channel.  A  threshold  is  determined  prior  to  proc¬ 
essing.  When  the  digital  count  in  channel  6  exceeds 
the  threshold,  its  location  is  determined  in  the  grid 
and  all  bits  are  set  denoting  land.  Data  points  with 
values  of  0  or  255  are  flagged  as  invalid.  Their 
position  on  the  grid  is  computed  and  the  high  bit  in 
the  word  is  set. 

At  the  completion  of  this  module,  each  pixel  can 
be  addressed  by  latitude  and  longitude  via  the 
mapping  equations. 

GEOMATCH.  Module  GEOMATCH  registers 
laser-  or  sonar-derived  control  points  to  the  rectified 
multispectral  image.  The  registration  is  accomplished 
by  using  the  mapping  coefficients  computed  in 
module  GEORECT. 

The  module  has  the  option  to  either  search  a 
portion  of  the  laser  sounder  file  or  search  the  entire 
file  for  soundings  that  lie  within  the  image  bound¬ 
aries.  If  laser  data  were  collected  concurrently  with 
the  multispectral  data,  then  only  that  portion  of  the 
file  is  searched.  A  time  interval  is  computed  based 
on  the  aircraft’s  speed,  altitude,  and  laser  angle  to 
insure  no  laser  sounder  depths  are  lost  due  to  the 
laser  sounder  leading  or  trailing  the  aircraft.  Since 
the  laser  sounder  file  is  sorted  by  time,  the  search 
time  is  significantly  reduced.  To  handle  the  case 
where  laser  sounder  data  may  have  been  collected 
at  different  times,  the  file  is  searched  sequentially. 


Other  data  sources  (i.e.,  data  collected  from  ships) 
can  also  be  used  by  this  module.  The  data  necessary 
are  latitude,  longitude  and  depth.  The  module 
matches  the  corresponding  pixel  values  based  on 
the  mapping  equations  and  outputs  a  data  set  for 
input  into  the  depth  modeling  programs. 

LMS.  LMS  uses  the  control  data  from  the 
GEOMATCH  program  to  compute  weights  for  depth 
estimation  (Martinez,  1988;  see  App.  B).  The  weights 
are  gridded  to  produce  a  weight  matrix,  which  is 
then  multiplied  by  a  pixel  matrix  to  produce  depths. 
In  the  case  of  more  than  one  channel,  several  ma¬ 
trices  are  computed.  The  gridding  of  the  weights  is 
a  lengthy  process  and  where  no  overlap  is  present 
on  the  image  edges,  gridding  of  the  weights  is  not 
recommended.  There  is  simply  not  enough  data  to 
accurately  grid  the  weights.  This  factor  is  a  problem 
with  the  LMS  algorithm. 

DEEP.  DEEP  uses  multiple  regression  methods 
(Clark,  1988;  see  App.  B)  to  compute  weights  for 
depth  estimation.  A  Cholesky  decomposition  is  used 
for  rapid  solution  of  a  system  of  simultaneous 
equations  and  to  produce  a  least-squares  solution 
for  the  depth  mapping  coefficients  These  coefficients 
are  multiplied  with  the  image  to  produce  depths. 

The  resulting  depths  are  output  in  image  file 
format  (described  in  App.  C).  The  depths  are 
multiplied  by  10,  so  that  0.1-m  precision  is  retained 
in  the  output.  The  depths  range  from  0  to  200, 
representing  depths  from  0  to  20  m.  The  present 
system  outputs  depths  to  a  1-m  grid,  although  user- 
requested  output  resolution  is  allowed. 

Results 

Data  Registration  and 
Rectification  Results 

Test  flights  with  GPS  coverage  were  based  upon 
a  Rockwell  Collins  GPS  1.023  MHz  clear  acquisi¬ 
tion  (C/A)  ranging  signal  code  receiver  provided 
by  the  Naval  Air  Development  Center  (NADC)  and 
installed  on  the  BIRDSEYE  P-3  aircraft.  The 
positioning  available  from  C/A  codes  is  less  accu¬ 
rate  than  the  10.23-MHz  precision  (P)  ranging  signal 
to  be  used  operationally.  This  reduction  in  aircraft 
positioning  accuracy  reduces  overall  registration 
accuracy  computed  from  test  data  sets. 

The  data  used  to  verify  the  attitude  correction 
algorithms  were  collected  over  Slidell,  Louisiana 
(near  New  Orleans).  Sidewalks  and  driveway 
intersections,  easily  located  on  the  multispectral 


Table  4.  Ground  truth  evaluation  of  corrected 
GPS  positions  from  aircraft  readings.  Results 
are  given  in  meters. 


Station 

X  Error 

Y  Error 

Total  Error 
(m) 

1 

18.4 

27.6 

33.1 

2 

25.7 

4.2 

26.0 

3 

26  4 

11.7 

28.9 

4 

24.2 

16.3 

29.2 

5 

13.2 

15.3 

20.2 

6 

11.7 

16.1 

19.9 

7 

2.2 

23.1 

23.2 

8 

13.2 

32.4 

35.0 

9 

22.6 

14.8 

27.0 

10 

21.3 

1.7 

21.4 

11 

13.9 

8.9 

16.5 

12 

0.7 

7.6 

7.6 

13 

6.6 

8.0 

10.4 

14 

2.2 

12.0 

12.2 

15 

6.6 

15.5 

16.8 

16 

14.7 

31.4 

34.7 

17 

7.3 

26.3 

27.3 

18 

7.3 

34.1 

34.9 

19 

19.1 

30.7 

36.1 

image,  were  chosen  as  ground  stations.  These  land¬ 
marks  were  well  distributed  throughout  the  image. 
Geographic  positions  at  the  ground  stations  were 
collected  using  a  Motorola  Mini  Ranger  Eagle  GPS 
receiver. 

The  pixel  location  of  the  ground  station  was  found 
on  the  multispectral  image.  The  corresponding 
latitude  and  longitude  were  computed  using  the  set 
of  mapping  coefficients  computed  for  the  image. 
Table  4  shows  the  results  at  19  locations.  The  X  and 
Y  errors  are  the  distances  in  meters  in  the  east- 
west  and  north-south  directions,  respectively, 
between  the  image-computed  position  and  the  ground 
station  GPS  reading.  The  total  error  is  the  straight 
line  distance  in  meters.  The  original  and  rectified 
images  are  shown  in  Figure  4. 

The  root-mean-square  (rms)  error  of  25.7  m  is 
based  on  the  deviation  of  the  computed  pixel  location 
from  the  averaged  GPS  location  of  that  pixel  on 
the  ground.  Many  sources  of  error  are  inherent  in 
this  value;  therefore,  this  rms  is  not  a  true 
measurement  of  the  model  error.  The  accuracy  of 
GPS  positions  depends  on  the  position  dilution 
of  precision  (PDOP).  The  PDOP  is  a  measure  that 
relates  the  configuration  geometry  of  the  GPS 
satellites  to  the  positioning  accuracy  (Wells,  1986). 
The  positioning  accuracy  is  a  product  of  the  PDOP 


and  the  measurement  accuracy.  Since  the  PDOP  was 
not  available  for  the  aircraft  data  and  was  not 
recorded  in  the  field  experiments,  the  actual 
positioning  accuracy  could  not  be  measured  for  this 
experiment.  The  measurement  accuracy  for  the  GPS 
ground  readings  was  measured,  however,  by 
recording  several  readings  while  the  receiver  was 
stationary  and  the  PDOP  was  not  changing.  An 
observation  was  that  the  high  variability  was  due 
to  differing  PDOPs  on  different  days  recorded  at 
the  same  location.  In  most  cases,  there  was  more 
day-to-day  variability  than  variability  due  to  the 
model  estimates. 

The  registration  algorithms  were  applied  to  data 
collected  over  shallow  waters  in  the  Florida  Keys. 
The  study  area  is  shown  in  the  dashed  box  in  Fig¬ 
ure  5.  The  laser  and  multispectral  image  registration 
in  Figure  6  is  for  a  500-m2  area  off  the  Eastern  Sambo 
Shoals.  Since  the  laser  positions  were  not  validated 
with  control  points,  the  accuracy  of  this  registration 
could  not  be  determined. 

Glint  Correction  Results 

Two  independent  data  sets  were  flown  on  differ¬ 
ent  days  off  Panama  City,  Florida,  and  Key  West, 
Florida.  Both  data  sets  were  collected  using  the 
NS.  Different  solar  elevation  angles  were  used  to 
study  glint  (Lingsch,  1989).  The  Panama  City  data 
were  collected  in  the  spring  of  1988,  and  the  Key 
West  data  were  collected  in  the  summer  of  1988. 
Figures  7  through  10  represent  true-color  images 
of  bands  2,  4,  and  5.  Figures  7  and  9  are  the  raw 
data.  Figures  8  and  10  are  the  corrected  data  for 
glint  removal.  In  each  case  the  deep-water  region 
used  for  the  model  was  a  subsampled  512x512  pixel 
scene  or  3500  data  points. 

The  Panama  City  corrected  data  set  (Fig.  8) 
exhibits  a  visual  improvement  over  the  raw  data. 
Several  pixels  on  the  far  right  were  not  corrected 
because  they  did  not  meet  the  criteria  for  a  glint- 
contaminated  pixel.  The  offshore  bar  was  preserved, 
whereas  much  of  the  “white  water,”  or  breaking 
waves,  was  corrected.  Some  of  the  vegetation  and 
roads  were  processed  as  water  pixels  because  the 
signal  in  the  near-infrared  images  (band  6)  over 
these  regions  fell  above  the  glint  and  below  the 
land  thresholds. 

The  Key  West  corrected  data  set  (Fig.  10)  reveals 
bottom  features  not  easily  apparent  in  Figure  9. 
These  features  appear  to  be  of  the  same  intensity 
as  similar  features  that  appear  in  the  nonglinted 
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Figure  4  Example  o]  SS  image  aver  Slidell,  before  (topi  and  after  (bottom)  geometric  rectification 


in 


Figure  5.  Study  area  south  of  Key  West,  represented  by  dashed  lines  on  1:40,000  scale  NOAA  chart  no.  11445.  Soundings 
are  shown  in  feet. 
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Figure  8.  Corrected  image  demonstrates  a  visual  improvement  The  breaking  waves 
were  also  corrected 


Figure  10  Bottom  feature  not  apparent  m  original  visible  image  band  after 
glint  correction 
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area.  This  is  an  indication  that  only  the  surface- 
reflected  signal  was  removed.  The  details  of  the 
glint  removal  are  shown  in  Appendix  B. 

Depth  Estimate  Results 

The  LMS  and  Regression  models  (App.  B)  were 
run  extensively  in  the  Key  West  area  of  Eastern 
Sambo  Shoals  (Fig.  6).  The  bathymetric  ground  truth 
was  collected  by  NAVOCEANO  in  April  1987.  The 
survey  track  line  spacing  was  25  m,  with  a  25-m 
interval  between  reported  depth  soundings. 
NAVOCEANO  provided  smooth  sheets  reduced  from 
the  survey  data  in  both  hard  copy  and  computer 
compatible  form. 

The  results  for  a  comparison  of  LMS  versus 
multiband  regression  methods  are  shown  in  Table  5. 
Although  the  LMS  did  relatively  well  in  areas  of 
constant  bottom  type,  it  did  not  do  as  well  with 
highly  varying  bottom  types.  One  reason  is  that  the 
nature  of  the  algorithm  makes  the  coefficients  adapt 
at  different  rates,  depending  on  the  model  param¬ 
eters  and  roughness  of  the  bottom.  When  the  bottom 
changes,  the  algorithm  adapts  the  coefficients  by 
the  amount  of  error  in  the  predicted  versus  the 
observed  depth.  These  large  error  spikes  contribute 
heavily  to  overall  error,  since  the  squared  error  is 
used  to  measure  rms  errors. 

Tables  6  through  8  give  results  for  the  multiband 
model  for  different  band  combinations  used  in  test¬ 
ing  the  model  estimates  against  the  NAVOCEANO 
ground  truth  data.  The  combination  of  bands  2  and 
4  gave  results  similar  to  the  five-band  model.  High 
correlation  in  the  visible  bands  is  the  main  cause, 
since  the  addition  of  extra  channels  does  not  add 
extra  information.  The  results  also  show  a  depth 
dependence  in  the  errors  (Fig.  11),  although  the 
percent  error  is  reasonably  constant.  Although 
the  error  in  shallow  water  was  close  to  0.3  m,  at  or 
below  3  m,  the  errors  beyond  3  m  were  above  the 
0.3-m  IHO  standard. 

The  boat  sounding  depths  were  used  as  control 
data  for  depth  estimation  from  the  multispectral  data; 
the  boat  depth  data  showed  better  correlation  with  the 
multispectral  depth  data  than  did  the  laser  sounder 
depth  data.  Also,  more  confidence  could  be  placed 
in  the  registration  of  the  MSS  data  with  the  boat 
data.  The  positional  accuracy  of  the  laser  data  was 
never  validated.  When  reliable  laser  data  becomes 
available  and  problems  are  resolved,  more  tests  will 
be  done.  Plots  of  the  MSS  data  versus  the  boat  and 


laser  sounder  data  are  compared  in  Figures  12  and 
13  respectively.  The  computed  depth  image  for  the 
Eastern  Sambo  Shoals  area  is  shown  in  Figure  14. 


CPU  Time 

Run-time  results  are  given  for  a  512  x  512  image 
processed  on  a  Digital  Equipment  Corporation 
VAX  8700  computer  at  NOARL  (Table  9).  This 
machine  has  a  processing  speed  of  approximately 
6  million  instructions  per  second  (mips). 


Summary 

A  software  package  that  computes  shallow- water 
depths  from  multispectral  data  combined  with  boat- 
or  laser-derived  soundings  was  developed.  The  best 
algorithms  did  not  produce  accuracies  within  the 
IHO  standards  of  0.3  m  at  depths  greater  than  3  m. 
The  horizontal  accuracies  depend  highly  on  the 
constellation  of  the  GPS  satellites,  since  high  PDOPs 
give  unreliable  estimates.  The  horizontal  accuracy 
derived  from  a  ground  control  test  gave  accuracies 
of  25  m,  but  the  PDOPs  were  too  high  to  give  a 
good  control  point  position  estimate. 

High  correlation  in  the  multispectral  wavebands 
causes  difficulty  in  separating  bottom-type  effects 
from  depth  effects  on  the  return  of  the  signal.  Since 
modeling  the  depths  from  the  multispectral  data 

Table  5.  Comparison  of  RMS  error 
for  multiband  algorithms  and  LMS 
results  using  NAVOCEANO  ground 
truth  data. 


LMS 

Multiband 

Control  set 

1  44  m 

1.01  m 

Test  set 

1.53  m 

1.07  m 

Table  6.  Multiband  model  results  using  NS  Bands 
1  through  5.  Errors  are  in  meters. 


Max 

Depth 

Regr. 

R* 

No. 

Callb. 

Points 

Test 

Mean 

Error 

RMS 

Error 

10 

0.701 

270 

269 

0.079 

1.270 

9 

0.692 

265 

265 

0.075 

1  252 

8 

0.670 

236 

234 

0.079 

1.143 

7 

0.602 

190 

190 

-0.0005 

0  986 

6 

0.633 

142 

140 

-0.097E-05 

0  855 

5 

0.603 

73 

73 

0.058 

0  691 

4 

0  631 

38 

40 

0.062 

0  544 

3 

0.881 

12 

12 

-0.954E-06 

0  323 
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Table  7.  Two  channel  results  from  multiband  model  using  combinations 
of  Bands  1,  2,  3  and  4.  Errors  are  in  meters. 


Depth 

Sensor/ 

Bands 

R2 

(Actual-Calculated) 
Mean  RMS 

No.  Calib. 
Points 

No.  Test 
Points 

0-3 

NS/ 13 

0  587 

-0.4167E-01 

0.320 

12 

12 

NS/ 14 

0.843 

-0.9537E-06 

0.322 

12 

12 

NS/23 

0.666 

-0.8333E-01 

0.311 

12 

12 

NS/24 

0.828 

-0.9537E-06 

0.3227 

12 

12 

0-4 

NS/13 

0.374 

-0.9537E-06 

0.5701 

38 

40 

NS/14 

0.605 

0.8750E-01 

0.5743 

38 

40 

NS/23 

0.170 

-0.1250E-01 

0.4998 

38 

40 

NS/24 

0.488 

0.5000E-01 

0.5679 

38 

40 

0-5 

NS/13 

0.457 

0.1027E-01 

0.7590 

73 

73 

NS/14 

0.593 

0.7877E-01 

0.6785 

73 

73 

NS/23 

0.264 

-0.1027E-01 

0.8516 

73 

73 

NS/24 

0.484 

0.7877E-01 

0.7225 

73 

73 

0-6 

NS/13 

0.460 

-0  6429E-01 

0.9790 

142 

140 

NS/14 

0.614 

-0.3929E-01 

0  8368 

142 

140 

NS/23 

0.398 

-0.2143E-01 

1.204 

142 

140 

NS/24 

0.604 

-0.1429E-01 

0.8833 

142 

140 

0-7 

NS/13 

0.430 

0.7894E-02 

1.168 

190 

190 

NS/14 

0  588 

-0. 1053E-01 

1.004 

190 

190 

NS/23 

0.384 

0  2631E-01 

1.220 

190 

190 

NS/24 

0.570 

0. 1579E-01 

0.9943 

190 

190 

0-8 

NS/13 

0  469 

0.3205E-01 

1.351 

236 

234 

NS/14 

0.656 

0  6624E-01 

1.186 

236 

234 

NS/23 

0.419 

0  1068E-01 

1.34 

236 

234 

NS/24 

0.633 

0.5555E-01 

1.15 

236 

234 

0-9 

NS/13 

0.502 

0  6509E-01 

1.459 

265 

265 

NS/14 

0.681 

0.688  7E-01 

1.293 

265 

265 

NS/23 

0.420 

0.4245E-01 

1.485 

65 

265 

NS/24 

0.638 

0.31 13E-01 

1.243 

265 

265 

0-10 

NS/13 

0.518 

0.6227E-01 

1.484 

270 

269 

NS/ 14 

0.690 

0.6970E-01 

1.301 

270 

269 

NS/23 

0.432 

0.4740E-01 

1.514 

270 

269 

NS/24 

0.646 

0.6041E-01 

1.269 

270 

269 

requires  as  many  bands  as  there  are  cover  types, 
more  independence  between  bands  or  less  cover 
types  is  necessary.  Since  the  independence  condi¬ 
tion  cannot  be  overcome,  a  reduction  of  cover  type 
is  more  desirable.  The  number  of  cover  types  can 
be  reduced  by  conducting  the  analysis  over  a  smaller 
area.  The  problem  with  using  a  smaller  area  is  that 
the  number  of  soundings  must  be  denser  to  provide 
adequate  information  to  the  model.  The  success  of 
the  system  is  more  probable  over  areas  with  little 
variability  in  bottom  cover. 

The  multispectral  data  are  of  considerable  value 
alone,  since  the  high  resolution  and  the  visibility 
through  the  water  column  provide  a  wealth  of 
information  about  the  nature  of  the  coastal 
environment.  The  success  of  the  glint-removal 


algorithms  is  valuable  as  a  tool  for  “seeing"  through 
highly  glinted  water  surfaces.  The  use  of  common 
classification  algorithms  is  useful  for  finding 
anomalies  in  these  types  of  data.  These  tools,  along 
with  the  GPS  positioning  data,  are  valuable  for 
updating  existing  hydrographic  charts. 

In  conclusion,  the  MIDAS  system  was  completed 
for  use  with  the  ABS  system.  The  multispectral 
imagery  can  now  be  processed  for  deriving 
bathymetry  in  clear,  shallow  waters  down  to  20  m. 
The  system’s  accuracy  varies  with  depth;  for 
example,  for  water  depth  less  than  3  m,  accuracy  is 
0.33  m,  and  at  a  water  depth  of  10  m,  accuracy 
is  1.6  m.  The  system  may  be  used  as  is  to  compute 
bathymetric  grids  from  multispectral  data  in 
conjunction  with  hydrographic  boat  survey  data  or 
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Table  8.  Two  band  results  using  combinations  of  Band  5  and  Bands 
1  through  4.  Errors  are  in  meters. 


Test  Depth  Residuals 

Depth 

Sensor/ 

Bands 

R2 

(Actual-Calculated) 
Mean  RMS 

No.  Calib. 
Points 

No.  Test 
Points 

0-3 

NS/35 

0.574 

-0.9537E-06 

0  3227 

12 

12 

NS/45 

0.756 

-0.4167E-01 

0.3200 

12 

12 

0-4 

NS/35 

0.165 

-0.9537E-06 

0  5244 

38 

40 

NS/45 

0,305 

-0. 1 250E-01 

0.5122 

38 

40 

0-5 

NS/35 

0.207 

-0.3767E-01 

0.7935 

73 

73 

NS/45 

0.179 

-0.2397E-01 

0.8111 

73 

73 

0-6 

NS/35 

0.211 

-0.6786E-01 

1.504 

142 

140 

NS/45 

0.311 

-0.1179E-00 

1.006 

142 

140 

0-7 

NS/35 

0.273 

-0  7896E-02 

1.217 

190 

190 

NS/45 

0.357 

-0  2632E-01 

1.153 

190 

190 

0-8 

NS/35 

0.322 

-0.1710E-01 

1.424 

236 

234 

NS/45 

0.430 

-0.2351E-01 

1.373 

236 

234 

0-9 

NS/35 

0.265 

0.3868E-01 

1.651 

265 

265 

NS/45 

0.362 

-0.9444E-03 

1.585 

265 

265 

0-10 

NS/35 

0.281 

0.1208E-01 

1.671 

270 

269 

NS/45 

0  348 

-0.4648E-02 

1  621 

270 

269 

Figure  II.  The  rms  error  of  depth  estimates  plotted  against 
depth  for  multiband  regression  model  used  with  boat 
soundings  and  five  bands  of  NS  data  in  Key  West  test  area. 


laser-derived  survey  data.  Boat  survey  data  were 
used  to  develop  this  particular  system,  since  laser- 
derived  depths  were  not  as  reliable. 

Recommendations 

In  the  future,  better  data  collection  systems  and 
algorithms,  along  with  the  complete  GPS  constel¬ 
lation,  will  become  available.  This  improved  data 
collection  ability  will  yield  measured  depths  closer 
to  IHO  standards,  for  depths  greater  than  3  m.  Also, 
data  with  less  noise,  better  horizontal  accuracy,  and 
better  ground  truth,  along  with  more  knowledge  about 
the  optical  parameters  of  the  water,  could  give  much 
better  bathymetry  estimates. 

Bathymetric  chans  were  derived  using  this  system 
for  the  Ocean  Venture  '90  exercise  over  Isla  De 
Vieques,  Puerto  Rico.  Charts  were  developed  from 
multispectral  imagery  to  support  amphibious  assault 
force  exercises  and  to  determine  the  usefulness  of 
such  products  in  amphibious  warfare.  There  are  many 
other  uses  for  this  type  of  product,  especially  if 
models  can  be  generalized  to  areas  for  which  no 
bathymetry  ground  truth  exists. 

Since  present  techniques  can  give  rms  errors 
between  11%  and  16  %  of  water  depth,  they  can  be 


17 
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Figure  12.  Plot  of  NS  band  ratio  (ch.  4  to  ch.2)  versus  depths  from  boat  sounding  data  in  Key  West  test  area. 
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Figure  13.  Plot  of  NS  band  ratio  (ch.  4  to  ch.  2)  versus  depths  from  laser  sounder  data  in  Key  West  test  area. 
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Figure  14  Depth  image  off  northwest  of  Eastern  Sambo  Shoals  area 
t lower  right  corner)  Color  bars  show  depth  contour  intervals  in  meters. 


Table  9  CPU  time  to  process  a 
512  x  512  image. 


Module 

CPU  Seconds 

Initialization 

00:00.33 

DEBLOCK.  PNAV 

00  1 1  30 

GEORECT 

00  58  41 

GEOMATCH 

00  02  84 

DEEP  (inc  output) 

00  18  43 

Total  Time 

01  32.41 

useful  in  situations  where  relative  depth  (c.g.,  the 
shallowest  versus  deepest  areas)  along  a  coastline 
needs  to  be  determined.  Such  products  can  also  be 
used  to  update  old  charts  with  information  on  new 
channels,  dredging,  sediment  deposition  areas, 
coastline  erosion,  harbors,  bridges,  and  construction. 
This  type  of  information  is  easily  available  at  a 
low  cost  from  multispectral  imagery  derived  depth 
and  enhanced  coastline  charts. 
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Appendix  A 


NOARL  Scanner  System  and  Specifications 
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tza  IHI/IO  BB  ON  SlOCSl 

Figure  A-l .  Schematic  diagram  ofNOARL  scanner. 
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Figure  A-2 .  Schematic  diagram  ofNOARL  scanner  optical  paths. 
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SPECIFICATIONS: 


9  CHANNELS 

-  6  VISIBLE 

-  3  IR 

SCAN  MIRROR 

-  ROTATING  OPTICAL  FLAT 

-  10.16  CM  DIAMETER 

-  ROTATION  RANGE: 

30  TO  165  RPS 

DIGITAL  ELECTRONICS 

-  3142  RESOLUTION  ELEMENTS 

PER  SCAN 

-  628  VIDEO  RESOLUTION 

ELEMENTS  PER  SCAN 

-  DIGITIZATION  LEVEL  8  BITS 


CALIBRATION  SOURCES 

-  2  BLACK  BODIES  AND  A  FIELD  FILLING 
INTEGRATING  SPHERE  FOR  EACH 
MIRROR  REVOLUTION 

-  SCANNER  DATA  WILL  BE  CALIBRATED  IN 

A  POSTPROCESSING  MODE  WITH 
LASER  SOUNOINGS 

SYSTEM  CHARACTERISTICS 

-  WEIGHT-275  LB 

-  POWER 

-  28  VDC 


Figure  A -3.  NOARL  scanner  specifications. 


Figure  A -4.  NOARL  scanner  operational  scenario. 
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Appendix  B 


Algorithms  and  Techniques 
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POSITIONING 

The  Global  Positioning  System  (GPS)  provides  latitude,  longitude,  and  heading  approximately  every 
second.  A  scan  rate  of  100  scan  lines  per  second  yields  a  positional  update  every  hundred  lines.  Positions 
are  linearly  interpolated  between  updates  as  a  function  of  time  to  produce  latitudes  and  longitudes  for 
individual  scanlines. 

GEOMETRIC  CORRECTIONS 

Corrections  for  roll,  pitch,  yaw  and  scan  angle  distortion  are  applied  to  the  GPS  derived  latitudes  and 
longitudes  to  derive  positioning  for  each  pixel. 

The  first  step  corrects  for  scan  angle  distortion.  As  the  distance  from  nadir  increases,  the  pixel 
footprint  size  increases.  The  distance  to  the  center  of  each  pixel  relative  to  nadir  is  computed  by 

Dx  =  TAN((N  *  IFOV)  +  IFOV/2)  *  Alt,  (1) 

where 

N  =  Pixel  number  (negative  left  of  nadir, 
positive  right  of  nadir); 

IFOV  =  Instantaneous  field  of  view  (0.002  Radians); 

Alt  =  Aircraft  altitude  (m); 

Dx  =  Distance  from  nadir  to  center  of  pixel  N  (m). 

If  the  scanner  is  not  toll  compensated,  the  roll  correction  is  applied  and  the  equation  becomes 

Dx  =  TAN((N  *  IFOV)  +  (IFOV/2)  -  ROLL)  *  Alt.  (2) 

A  positive  roll  angle,  or  right  bank  roll,  translates  the  scanline  to  the  left,  and  a  negative  roll  translates 
the  line  to  the  right. 

The  pitch  correction,  Dy,  translates  the  entire  line  forward  or  backward  by  the  ground  distance  traveled 
when  the  scanner  rotates  about  the  angle  of  pitch. 

Dy  =  Alt  *  TAN(PITCH)  (3) 

The  offsets  are  rotated  by  the  true  heading  relative  to  the  geographic  coordinate  system.  The  true 
heading  is  computed  from  the  heading  received  from  the  GPS  and  the  drift  or  yaw  received  from  the  INS. 
The  GPS  heading  is  computed  from  point  to  point  locations.  The  drift  angle  is  then  defined  as  the 
deviation  of  the  aircraft’s  longitudinal  axis  from  the  GPS  heading.  A  positive  drift  means  the  nose  of  the 
aircraft  points  left  The  true  heading  is  given  by: 

8  =  Heading  -  Drift  (4) 

Each  pixel  is  then  multiplied  by  a  rotation  matrix,  equation  5,  and,  offsets  Dx’  and  Dy’  in  meters  are 
computed  for  the  pixel’s  distance  from  the  center  position. 

Dx]  [cose  SIN0  Idx’ 

.  =  (5) 

Dy  -SIN©  COSe|  [Dy’ 
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The  latitude  and  longitude  are  computed  using  conversion  formulas  from  Bowditch  (1975). 

Degrees  /  Meter  of  Latitude 
=  111132.09  -  566.05  *  COS(2  *  Latitude)  + 

1.2  *  COS(4  *  Latitude)  -  0.002  *  COS(6  *  Latitude) 

Degrees  /  Meter  of  Longitude 

=111415.13  *  COS(Latitude)  -  94.55  *  COS(3  *  Latitude)  + 

0.12  *  COS(5  *  Latitude) 

SUBAREA  DETERMINATION  AND  MAPPING  EQUATIONS 

Each  continuous  flight  line  is  divided  into  subareas  for  processing.  Each  subarea  covers  approximately 
a  600  x  500  m  area.  The  first  step  computes  the  boundaries  of  the  subarea.  Distances  AB  and  DE  in 
Figure  B-l  are  computed  by 

Alt  *  TAN((Pix/2  -  Ste)  *  IFOV),  (6) 


where 


Pix  =  Number  of  pixels  in  scanline. 

Ste  =  Starting  element  number. 

Distances  AC  and  DF  are  computed  by 

Alt  *  TAN((Ene  -  Pix/2)  *  IFOV),  (7) 

where 

Ene  =  ending  element  number. 

Variables  Ste  and  Ene  are  region-of-interest  boundaries,  which  define  the  area  to  be  processed. 
Portions  of  the  image  may  not  be  processed  due  to  sun  glint/glare  or  other  environmental  factors. 

Once  these  distances  are  computed,  B  and  C  are  rotated  about  A,  and  E  and  F  are  rotated  about  D  by 
the  heading.  B\  C\  E’,  and  F’  become  the  subarea  comer  points.  Approximately  512  scanlines  are 
processed  in  a  subarea;  hence,  D  is  about  512  lines  from  A.  This  size  depends  on  the  scan  rate  and  the 
aircraft’s  speed.  A  and  D  are  known  geographic  locations  that  have  been  corrected  for  aircraft  attitude. 
Geographic  locations  are  computed  for  the  comer  points.  The  top  boundary  of  the  subarea  will  become 
the  lower  boundary  of  the  subsequent  subarea  processed. 

Six  control  points  are  defined  as  shown  in  Figure  B-l  (A,  B’,  C\  D,  E\  and  F’).  Additional  control 
points  are  added  by  sampling  every  50th  scanline.  The  distance  of  the  center  pixel  from  the  center  of  the 
subarea  is  computed  by  using  the  Law  of  Cosines.  Twenty  control  points  along  the  scanline  are 
geometrically  corrected  as  in  the  previous  step,  giving  roughly  200  equally  spaced  control  points  within 
the  subarea.  Once  all  the  control  points  are  computed,  mapping  coefficients  are  computed  using  a  second 
order  polynomial.  The  following  equation  converts  latitude  and  longitude  to  line  and  element: 

Line  =  B,  +  B,  Lat  +  B3  Lon  +  B4  Lat2  +  B,  Lon2  +  B6  Lat  Lon,  (8) 

Elem  =  C,  +  Q  Lat  +  C3  Lon  +  C4  Lat2  +  Cs  Lon2  +  C4  Lat  Lon,  (9) 
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where 


B„  =  Line  coefficients, 

Cn  =  Element  coefficients. 


E 


Figure  B-l.  Figure  shows  rotation  of  subvea  for  heading.  Lines  BC  and  EF  are  rotated,  B',  C\ 
and  F‘  become  comer  boundaries  of  subarea. 

Each  pixel  goes  through  the  transformation  for  aircraft  attitude  and  is  mapped  to  the  geographic  grid 
by  equations  8  and  9.  Data  gaps  occur  for  the  following  reasons.  As  the  distance  from  nadir  increases, 
the  distance  between  the  center  of  pixels  increases,  i.e.,  at  an  altitude  of  500  m  the  distance  between  pixels 
becomes  1.7  m  at  the  farthest  off-nadir  positioa  When  transformed  to  a  1-m  grid,  this  effect  causes  a  gap 
on  the  vertical  axis.  Gaps  in  the  horizontal  axis  occur  depending  on  the  sample  rate,  air  speed,  and 
changes  in  the  aircraft’s  attitude.  Nearest  neighbor  interpolation  is  used  to  fill  in  these  gaps.  If  the  vacant 
cell  does  not  have  an  adjacent  neighbor,  no  value  is  placed  in  that  cell  and  no  depth  is  computed  for  that 
position. 

Equations  8  and  9  can  be  used  to  register  any  external  source  of  soundings  to  the  image.  An 
advantage  to  this  method  is  any  data  set  can  be  registered  to  the  multi  spectral  image  regardless  of  the 
source  or  the  time  collected.  At  die  same  time  die  coefficients  are  generated,  coefficients  for  the  inverse 
problem  are  computed  (i.e.,  given  line  and  element  compute  latitude  and  longitude).  This  is  useful  for 
later  reference  and  mosiacking  subareas  to  produce  a  complete  bathymetric  chart. 
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GLINT  REMOVAL 

Passive  multispectral  data  collected  over  water  are  subject  to  both  specular  and  diffuse  reflection  off 
the  surface  of  the  water.  This  component  of  the  signal  appears  as  glint  or  glare  at  the  sea  surface  and 
must  be  removed  from  the  signal  for  bathymetry  estimators  to  be  effective. 

Several  authors  (Lyzenga,  1985;  Spitzer,  1987)  have  reported  the  signal  in  the  near  infrared  (IR)  band 
is  useful  in  determining  the  magnitude  of  surface  reflection.  Transmittance  of  visible  wavelengths  is  good, 
but  a  significant  absorption  band  occurs  in  the  near-IR  band  (750-760  nm)  (Jerlov,  1968).  As  a  result, 
while  surface  reflection  and  atmospheric  returns  are  common  to  all  bands,  the  light  reflected  from  the 
water  column  and  the  bottom  is  mostly  light  in  the  visible  wavelengths  (450-690  nm). 

Lyzenga  (1985)  proposes  a  method  that  exploits  the  spectral  characteristics  of  sun  glint  to  derive  a 
correction  for  surface  reflection.  He  assumes  a  linear  relationship  between  the  glint  in  the  visible  and  IR 
channels.  This  results  in  a  glint  correction  for  visible  band  i  based  on  IR  channel  j  of  the  form, 

AYj  =A(Yj-Yj).  (10) 

The  corrected  visible  band  signal  is  given  by 

%  =Yi-AYi,  (11) 

where 

Y,  =  Raw  signal  for  visible  band, 

Yj  =  Average  IR  signal  in  a  non-glinted  deep  water  region, 

oy  =  Covariance  of  the  visible  and  the  near  IR  bands 
computed  in  a  deep  water  region, 

Oy  =  Variance  of  the  near-IR  band, 

A  =  Oy  /  Oy. 

Lyzenga  suggests  the  method  is  valid  when  the  correlation  coefficient  is  greater  than  0.90. 

Two  observations  motivate  modifications  and  extensions  of  Lyzenga’s  method.  A  linear  model  may 
not  accurately  represent  the  relationship  between  the  output  of  the  scanner  sensors.  Figure  B-2  is  a  typical 
scatter  plot  of  raw  data  from  the  NOARL  scanner  channel  6  versus  channel  4.  Motivation  for  a  higher 
order  model  is  clearly  given  in  the  nonlinear  relationship  between  the  channels.  Second,  in  shallow  water, 
the  indiscriminate  application  of  Lyzenga’s  method  has  been  seen  to  remove  prominent  features  such  as 
sandbars.  The  solution  proposed  here  is  to  identify  and  mask  the  glint-free  areas  in  the  image  and  to 
process  only  those  areas  with  glint 
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DEEP  WATER  STATISTICS 

The  procedure  for  glint  removal  first  involves  preprocessing  a  deep  water  region  containing  both  glint 
and  nonglint.  A  deep-water  region  is  chosen  to  assure  that  the  signal  in  the  visible  bands  is  independent 
of  bottom  reflection.  Since  the  IR  does  not  penetrate  the  water  column,  the  covariance  in  the  IR  and 
visible  channels  over  deep-water  is  assumed  to  be  due  to  glint. 

Glint-contaminated  regions  exhibit  a  large  pixel-to-pixel  variance,  and  nonglint  regions,  a  smaller 
variance.  This  variability  is  clearly  seen  in  Figures  B-3  and  B-4,  which  show  a  typical  signal  cross  section 
of  a  deep-water  region  in  the  near  IR  band  and  the  corresponding  image.  Using  this  property,  a  nonglint 
area  within  the  deep-water  region  is  found  and  the  averages  Yj  and  Y*  are  computed. 

A  mask  is  produced  to  insure  that  only  areas  with  glint  are  processed  and  that  such  features  as  land, 
nonglint  areas,  and  sandbars  are  not  altered. 

The  near-IR  band  6  is  used  to  distinguish  land  and  water.  A  threshold  is  selected  manually.  An  area 
in  the  flight  line  is  chosen  that  contains  both  land  and  water.  By  highlighting  pixels  in  the  near  IR  from 
0  to  255  successively,  a  maximum  water  value  is  obtained.  Ideally,  IR  band  7  could  be  expected  to 
distinguish  surface-subsurface  reflections  better  than  band  6.  However,  band  7  was  excessively  noisy  and 
failed  to  sufficiently  resolve  the  land-water  boundary. 

The  IR  channel  is  also  used  to  distinguish  between  possible  glint  and  nonglint  pixels.  Over  water,  the 
near-IR  band  is  assumed  to  consist  almost  entirely  of  surface  reflection.  Thus,  a  strong  return  is  an 
indication  of  glint  in  a  pixel. 


Figure  B-3.  Channel  6  cross  section  for  line  A-B  in  Figure  5  below,  showing  the  transition  from  a 
nong tinted  region  to  a  glinted  region. 
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Figure  B-4.NS  near-IR  band (700-900  run).  Signal  over  shoals  and  sandbars  showing  tone,  whereas  the  signal 
due  to  reflection  demonstrates  texture. 


QUADRATIC  GLINT  CORRECTION 

The  method  discussed  here  is  a  variation  of  Lyzenga’s  method,  with  the  inclusion  of  a  squared  term 
in  the  model.  The  glint  correction  for  the  visible-band  i  based  on  the  IR-band  j  takes  the  form, 

AY  =  A(Yj  -  Y;)  +  B(Y2;  -  Y2 ).  (12) 

i 

The  corrected  visible-band  signal  is  given  by 

Y,  =  Y,  -  AY,,  (13) 

where 

Y,  =  Raw  signal  for  visible  band, 

Yj  =  Raw  signal  for  ncar-IR  band, 

Yj  =  Average  IR  signal  in  a  nonglint  deep-water  region, 

Y2  =  Average  IR  signal  squared  in  a  nonglinl  deep-water  region,  and 

A,  B  =  Regression  coefficients. 

Application  of  this  model  requires  the  identification  of  suitable  non-glint  deep  water  regions  over  which 
the  averages  Y;  and  Y2  arc  computed.  To  compute  the  regression  coefficients,  each  of  the  five  visible 
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bands  is  regressed  against  the  near-IR  band.  The  region  where  the  coefficients  are  computed  requires  a 
deep-water  region  containing  both  glint  and  nonglinL 

DEPTH  ESTIMATION  REGRESSION  MODEL  ESTIMATION 

Current  depth  algorithms  are  based  on  the  functional  form  of  the  single-band  radiance  model  (Clark 
et  al.,1987).  Statistical  regression  techniques  are  used  to  determine  the  unknown  parameters  using  water 
depths  known  at  a  few  pixel  locations  in  die  multispectral  image.  Single-  and  dual-band  ratio  models  are 
very  sensitive  to  changes  in  bottom  reflectance  and  water  clarity.  Multiband  algorithms  in  which  the 
radiance  and  bottom  reflectance  are  functions  of  the  wavelength  can  be  used  to  minimize  difficulties 
encountered  due  to  changes  in  the  bottom  reflectance  and  water  attenuation.  A  development  of  the 
multiband  model  from  the  single-band  radiance  equation  is  presented  here. 

The  single-band  radiance  equation,  where  total  radiance  from  the  water,  Ljt  for  a  single  bandwidth  (i), 
is  written  as 

Li  =  Lmi  +  Si  Tj  exp[-f  lq  z],  (14) 


where 


L_  =  radiance  from  deep  water, 

S;  =  a  function  of  transmittance  of  the  air  and  water  surface,  surface  refraction  and  solar 
irradiance; 

r(  =  bottom  reflectance; 

=  effective  attenuation  coefficient  of  the  water, 
f  =  geometric  factor  to  account  for  the  pathlength  through  the  water, 
z  =  water  depth; 
i  =  l,2...n; 

n  =  number  of  spectral  bandwidths. 

There  have  been  several  implementations  of  this  model  (Gark  et  al.,1987;  Paredes  and  Spero,1983).  Some 
implementations  require  classifying  the  data  into  bottom  cover  types  or  knowledge  of  optical  parameters. 
Other  implementations  use  band  ratios  and  combinations  of  bands.  Ratios  work  rather  well  if  they  are 
constant  for  all  bottom  types  in  the  area  of  interest.  The  ratio  method  was  outlined  by  Polcyn  et  al.(1976) 
and  elaborated  by  Lyzenga  (1985).  The  generalized  band  ratio  algorithm  has  been  found  to  give  the  best 
results  as  discussed  by  Gaik  et  al.  (1988).  A  brief  summary  of  the  method  is  given  here.  Equation  (14) 
is  usually  linearized  and  solved  for  z  to  yield 

z  =  (1/fk,)  (lntrj  -  X,).  (15) 

where  X,  =  ln((L,  -  L^/sJ.  Equation  (2)  can  be  reparameterized  with 


b*  =  (1/flq)  (ln[rj), 
bn  =  -(1/flq) 
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and 

X,,=  l  . 

The  depth,  z,  for  one  channel,  can  then  be  expressed  in  matrix  notation  as 
z  =  b’x, 

where  b’  is  the  transpose  of  the  vector  b.  The  vector  b  contains  the  weights  bo  and  bj.  The  vector  x 
contains  [1.x], 

For  n  bandwidths,  the  depth  z  can  be  expressed  as  a  weighted  sum  of  single-band  terms  as  follows 
z  =  w,  b,  x,  +  w2  b*  Xj  +  ...  +  w„  bB  x„, 

where  the  b  and  x  vectors  are  defined  as  above  for  each  band,  and  the  w’s  are  weights  whose  sum  equals 
1.  In  matrix  form,  this  equation  can  be  written  as 

z  =  w’B’x,  (16) 

where  w  is  a  vector  of  n+1  weights,  B  is  a  diagonal  matrix  of  weights  where  bo  is  a  constant  term 
derived  from  the  stun  of  constant  terms  for  each  band,  and  x  is  a  n+1  vector  with  Xo=  1.  The  motivation 
behind  the  summation  of  constant  terms  (Paredes  and  Spero,  1983),  is  discussed  here. 

For  m  bottom  cover  types,  and  n  bandwidths  such  that  m  =  n,  Paredes  and  Spero  have  shown  that  a 
vector  c  of  length  n  exists,  such  that 

c’A  =  1, 

where  1  is  a  l’s  vector  of  length  n  and  A  is  defined  as  the  n  by  m  matrix. 


Each  element  of  A  is  the  logarithm  of  the  bottom  reflectance  term  for  the  nth  channel  and  the  mth  bottom 
type.  As  long  as  the  number  of  bottom  types,  m,  is  not  greater  than  the  number  of  channels,  n,  a  solution 
for  c  can  be  found.  Spero  and  Paredes  show  that  for  m=n  bottom  types,  the  weights  wt  can  be  written 

w,  =  Cilq  /  c’k  . 

Equation  16,  can  be  written  as 
z  =  who  +  W’M, 

where 

w’bo  =  (l/Otc’kr1  *  c’A. 
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Since  c’A  is  1,  the  bottom  reflectance  term  goes  away,  and  only  the  sum  of  constant  terms  remains.  The 
multiband  depth  model  can  thus  be  expressed,  independently  of  bottom  type  reflectance,  as  the  linear 
regression  model 

z  =  px  +  e, 

where  0  is  the  vector  w’B  of  coefficients  derived  from  the  least-squares  solution  of 
equation  16,  X  is  the  vector  of  reflectances,  and  e  is  the  error. 

Results  from  comparison  of  the  multiband  model  to  the  single-band  and  dual-band  ratio  modeL 
indicate  that  the  multiband  generalized  ratio  technique  gives  superior  results  and  is  computationally 
efficient  as  well. 

THE  LMS  ALGORITHM 

Another  approach  taken  is  an  adaptive  fitting  technique,  often  used  in  linear  systems  theory,  called 
the  LMS  algorithm  (Widrow,  1985). 

Depths  are  estimated  by 

i  =  w’x, 

where  i  =  estimated  depth; 

w  =  vector  of  weights,  (w’  is  transpose); 
x  =  vector  of  reflectance  values. 

The  error  at  each  point  is  given  by 
e  =  z  -t. 

The  weights  are  computed  at  each  known  depth,  as  follows: 
wk+,  =  wk  +  2\x&^, 


where 


wk  =  weight  for  k*  point; 

p.  =  constant  chosen  for  stability  and 
rate  of  convergence; 

ej.  =  error  for  the  k*  point; 

xk  =  input  value  for  the  k®  point. 

For  more  than  one  bandwidth,  the  equation  may  be  expressed  in  matrix  notation, 
wk+,  =  wk  +  2pex, 
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where  wk  is  a  vector  of  weights  for  each  band,  and  x  is  a  vector  of  reflectances  in  each  band.  Optionally, 
the  x’s  may  be  ratios,  and  a  vector  of  ones  can  also  be  included  to  add  stability. 

Computing  the  weights  is  quite  simple  and  does  not  require  much  computation  time.  Preliminary 
results  have  shown  that  this  technique  can  produce  more  accurate  results  than  the  use  of  a  deterministic 
model. 
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Appendix  C 


Depth  File  Format 
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DEPTH  DATA  FILE  FORMAT 

A.  File  Name 

DISK$SCANNER:[MIDAS.DATA.Mnn]TnnFnn.DEP  where: 

Mnn  =  The  survey  mission  number  (ie.,  M12) 

Tnn  =  The  tape  number  (ie.,  T24) 

Fnn  =  The  file  number  (ie.,  F02) 

B.  Usage 

This  file  is  the  end  product  of  the  Multispectral  Image  Depth  Analysis  System  (MIDAS).  This  file 
contains  all  the  necessary  information  needed  to  generate  a  navigation  chart  for  a  specified  survey  area. 

C.  File  Structure 

The  DEPTH  data  file  is  a  variable-length,  random-access,  block-structured  file.  It  contains  a  main  header, 
which  describes  general  information  about  the  data  an  survey  area.  All  depth  data  are  stored  via  linked 
list  subarea  headers,  which  contain  information  about  the  subarea  and  a  pointer  to  the  first  block  of  the 
packed  depth  data  for  the  subarea.  Refer  to  the  next  page  for  a  description  of  a  typical  MIDAS  Depth 
Data  File  format. 

D.  Record  Structure 

There  are  three  primary  record  types  in  the  MIDAS  Depth  Data  File:  The  Main  Header,  The  Sub-Area 
Headers,  and  the  Packed  Depth  Data.  Refer  to  the  following  pages  for  detailed  descriptions  of  each  record 
structure. 


Record  fl 

Pointer  to  1st  j 
Subtree  Header.  ' -> 
.-> 

The  subarea 

header  linked 
list  pointers 
(Next, (.Last) 
link  the 

subarea  headers 

together. 


Linkage  between 
the  2nd  SubArea 
Header  and  the 
3rd  SubArea 
Header. 


Link  between  j 
SubArea  Headers  v 
3  and  4. 


Main 

Header 


1st  SubArea 

Header 


”i 


Packed  Depths  For 
The  1st  SubArea 


<-> 


2nd  SubArea 

Header 


Packed  Depths  For 
The  2nd  SubArea 


3rd  SubArea 

Header 


<-• 


Pointer  to  the 
first  block  of 
packed  depths 
for  SubArea  1. 


Pointer  to  the 
first  block  of 
packed  depths 
for  SubArea  2. 


Pointer  to  the 
first  block  of 
packed  depths 
for  SubArea  3. 


Itc... 


!-> 


.  -> 


Note  that  "Next" 
in  the  last 
SubArea  Header 
points  to  the 
next  free  block 
in  the  file. 

The  nuaber  of 
SubArees  in  the 
file  is  found 
in  the  Main 
Header. 


-> 


Last  SubArea 
Header 


Packed  Depths  For 
The  Last  SubArea 


Next  Free  Slock 
Data  File 


In 


Pointer  to  the 
first  block  of 
packed  depths 
for  the  last 
SubArea  in  the 
file. 
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E.  Main  Header 

Record  Type:  Main  Header  (For  MIDAS  Version  #1) 

Record  Access:  Block  I/O  (suggest  ELTRAN) 

Record  Structure:  Fixed  length  formatted  (1024)  byte)  block. 

Of  the  1024  total  bytes,  409  bytes  are  used  with  615  bytes  reserved  for  expansion.  The  Version  Number 
(Vers)  of  the  file  at  offset  4  defines  what  record  structures  to  use  when  reading  the  depth  file.  At  the  time 
this  documentation  was  written  there  was  only  one  format  (Vers  1  format).  This  variable  should  be  used 
by  the  software  to  keep  the  software  upward  compatible.  The  contents  of  the  Main  Header  block(s)  are 
as  follows: 


Offset 

(bytes) 

Name 

Size 

(bytes) 

Date 

Type 

Description 

0 

Idem 

4 

CHAR 

File  type  identification  code  ‘DEEP' 

4 

Ven 

4 

LONG 

Version  number  of  creating  software 

8 

SunAOf 

4 

LONG 

Offset  to  first  subarea 

12 

NoSubA 

4 

LONG 

Number  of  subareas  in  the  file 

16 

Mission 

4 

LONG 

Mission  number 

20 

Tape 

4 

LONG 

Tape  number 

24 

File 

4 

LONG 

Number  of  file  on  the  tape 

28 

Resoln 

4 

REAL 

Resolution  of  the  daU  (in  meters) 

32 

MinLat 

8 

DBLE 

Minimum  latitude  value  in  this  file 

40 

MinLon 

8 

DBLE 

Minimum  longitude  value  in  this  file 

48 

MaxLac 

8 

DBLE 

Maxim  tan  latitude  value  in  this  file 

56 

MaxLon 

8 

DBLE 

Maximan  longitude  value  in  this  file 

64 

Date 

9 

CHAR 

Date  of  survey 

73 

Area 

80 

CHAR 

Name  of  the  area  coveted 

153 

Content 

255 

CHAR 

Additional  comments  (of  any  type) 

408 

IsQuad 

1 

FLAG 

Flag  indicating  if  quadratic 
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Data  Type 
CHAR 


Size  (bytes) 


VAX  Fottnn  77  Equivalent 
Character  *  4 


CHAR 

80 

Character  *  80 

DBLE 

8 

Real  •  8 

DBLE 

32 

Real  *  8  array  of  4 

FLAG 

1 

Logical  *  1  (.true,  or  iiise.) 

LONG 

4 

Integer  *  4 

REAL 

4 

Real  *4 

F.  Sub  Area  Header 

Record  Type:  SubArea  Header  (for  MIDAS  Version  #1) 

Record  Access:  Block  I/O  (suggest  ELTRAN) 

Record  Structure:  Fixed  length  formatted  (512  byte)  block. 

Of  the  512  total  bytes,  192  bytes  are  used  with  320  bytes  reserved  for  expansion.  The  contents 
of  the  Subarea  Header  Block  are  as  follows: 


Offset 

(bytes) 

Name 

Size 

(bytes) 

Date 

Type 

Description 

0 

Next 

4 

LONG 

Block  number  of  next  subarea  header 

4 

Last 

4 

LONG 

Block  number  of  last  subarea  header 

8 

DephOf 

4 

LONG 

Starting  block  number  of  depth  data 

12 

Nun)  Sen 

4 

LONG 

Number  of  scanNumSea  in  subarea 

16 

NumHm 

4 

LONG 

Number  of  elements  in  subarea 

20 

BgnCov 

4 

LONG 

Beginning  HALS  coverage  element 

24 

EndCov 

4 

LONG 

Ending  HALS  coverage  element 

28 

SbALat 

32 

DBLE 

Comer  Latitudes  for  subarea 

60 

SbAlon 

32 

DBLE 

Comer  longitude  for  subarea 

92 

UnCof 

48 

DBLE 

The  scanline  mapping  coefficients 

140 

EhnCof 

48 

DBLE 

The  elements  mapping  coefficients 

188 

SbArea 

4 

LONG 

The  physical  Subarea  number 
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Data  Type 
CHAR 


Size  (bytea) 


VAX  FORTRAN  77  Equivalent 
Chancier  *  4 


CHAR 

80 

Chancier  •  80 

DBLE 

8 

Real*  8 

DBLE 

32 

Real  *  8  array  of  4 

FLAG 

1 

Logical  *  1  (.true,  or  .false.) 

LONG 

4 

Integer  *  4 

REAL 

4 

Real*  4 

G.  Packed  Depth  Format 

Record  Type:  Packed  Depth  Data  (for  MIDAS  Version  #1) 

Record  Access:  Block  I/O  (suggest  MIDAS  Pack/Unpack  Routines) 

Record  Structure:  Scaled  decimal  (byte)  data  packed  into  a  contiguous  set  of  (512  byte 

blocks. 

The  last  block  in  the  set  is  padded  with  ASCII  nulls. 


Example  of  511  byte  records  packed  into  512  byte  blocks: 


BYTE# 

BYTE# 

H 

BYTE# 

BYTE# 

BYTE# 

BLOCK 

1 

2 

H 

510 

511 

512 

1 

1.  1 

2.  1 

H 

510,  1 

511,  1 

1.  2 

2 

2.  2 

3.  2 

511,  2 

1.  3 

2.  3 

3 

3.  3 

4,  3 

■ 

1.  4 

2.  4 

3.  4 

•  a. 

... 

... 

... 

... 

510 

510,  510 

511, 510 

■ 

509,511 

510,511 

511 

511,  511 

1,512 

mmm 

510,  512 

511,  512 

512 

1,513 

2,513 

B 

510,  513 

511,  513 

1,514 

513 

2,514 

3,514 

B 

511,  514 

<BLANK> 

<BLANK> 

ELEMENT.  SCANUNE 

Each  Mock  represents  (depth  X  10)  stored  la  a  stag! e  byte  at  the  specified 
Elen  eat  and  ScaaHae. 

This  example  shows  the  typical  storage  of  the  packed  depth  values  for  a  511  element  by  514 
scanline  subarea.  There  are  no  extra  characters  (such  as  <CRxLF>)  stored  in  the  packed  data. 
The  pack  and  UnPack  routines  use  the  number  of  elements  and  number  of  scanlines  stored  in  the 
subarea  header  to  pack  and  unpack  the  data.  Depth  time  10  values  are  stored  in  bytes.  Valid 
depths  are  0.1  -  20.0  meters  stored  as  1-200  in  the  file.  A  zero  (0)  represents  a  land  value,  and 
a  255  represents  an  invalid  data  point.  The  last  block  written  will  be  added  with  binary  zero  (0) 
or  ASCII  null. 
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